function runAnalysis(Params)

% This is the main function that performs ISC analysis for 
% preprocessed fMRI data. Analysis parameters (see initParams.m) 
% must be set by the user before running the analysis.
%
% Input fMRI data must be either in "nii" or "mat" -format 
% (see initParams.m). Output data is mapped into the memory and 
% can be quickly accessed through memory pointer objects.
%
% COMMAND LINE ANALYSIS:
% 
% After running the analysis, memory pointer objects can be loaded 
% to the Matlab's workspace by typing:
% load memMaps % objects are saved in analysis destination folder.
%
% GRAPHICAL USER INTERFACE (GUI):
% 
% Specific GUI has been designed to visualize and analyze the results.
% After running the analysis, GUI can be launched by typing:
% load Tag % Tag is saved in the analysis destination folder
% ISCtool(Tag); % launch GUI from your GUI folder
% 
% Note: when you launch GUI make sure you have cleared all memory
% map pointer objects from Matlab workspace.
% 
% See also: ISCANALYSIS, MEMMAPDATA, INITPARAMS

% Last updated: 2.9.2010 by Jukka-Pekka Kauppi
% Tampere University of Technology
% Department of Signal Processing
% e-mail: jukka-pekka.kauppi@tut.fi


% BATCH 1:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Initialize parameters. Note: user must set all public parameters
% in initParams.m before running the analysis.
%Params = initParams;
% initialize data matrices into the disk and create pointers
% using dynamic memory mapping:
Priv = Params.PrivateParams;
Pub = Params.PublicParams;

Params = memMapData(Params);

% BATCH 2:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Filter all data through stationary wavelet filter bank.
% Computationally faster way is to distribute the
% calculations across several processors.
disp(' ')
disp('Filtering:')
for nrSubject = 1:Priv.nrSubjects
  for nrSession = 1:Priv.nrSessions
    disp(['Subject: ' num2str(nrSubject) ...
    ' , Session: ' num2str(nrSession) ':'])
    filterData(Params,nrSubject,nrSession);
  end
end

% BATCH 3:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Calculate intersubject synchronization maps for every frequency band.
% Computationally faster way is to distribute the
% calculations over several processors.
% Here, nrBand 0 stands for full frequency band.
disp(' ')
disp(['Calculating inter-subject synchronization maps:'])
for nrBand = 0:Pub.nrFreqBands
  for nrSession = 1:Priv.nrSessions
    disp(['Band: ' num2str(nrBand) ', Session: ' num2str(nrSession)])
    calculateSimilarityMaps(Params,nrBand,nrSession);
  end
end

% Assess statistical significance within each frequency-specific
% ISC map through permutation testing.
% Generate permutation (null) distributions:
disp(' ')
disp('Calculating permutation distributions for each ISC map:')
for nrSession = 1:Priv.nrSessions
  for permSetIdx = 1:Priv.nrPermutationSets
    permutationTest(Params,nrSession,permSetIdx,0); % across session
    permutationTest(Params,nrSession,permSetIdx,1); % time-windows
  end
end

% Calculate frequency comparison maps (sum ZPF statistic).
% Get total number of frequency band comparisons:
disp(' ')
disp('Calculating sum ZPF statistic for frequency-band comparisons:')
freqComps = ((Priv.maxScale+2)^2-(Priv.maxScale+2))/2;
%  Calculate null distributions for all comparisons indexed from
% 1 to freqComps. To distribute computations, you can also call
% function  by giving subblocks as input, e.g. call function
% separately with comparisons 1:3, 4:6, 7:9, 10:12, and 13:15.
for nrSession = 1:Priv.nrSessions
  PearsonFilon(Params,nrSession,1:freqComps);
end

% Assess statistical significance between frequency bands through
% permutation testing.
disp(' ')
disp('Calculating permutation distributions for sum ZPF maps:')
for nrSession = 1:Priv.nrSessions
  for freqComp = 1:freqComps
    permutationPF(Params,nrSession,freqComp)
  end
end

% Calculate full intersubject correlation matrices. These calculations
% are required if subject-wise correlations are investigated 
% (average ISC is calculated using the function 
% calculateSimilarityMaps.m but it saves only the mean correlation values).
disp(' ')
disp('Calculating inter-subject correlation matrices:')
for nrBand = 0:Pub.nrFreqBands
  for nrSession = 1:Priv.nrSessions
    disp(['Band: ' num2str(nrBand) ', Session: ' num2str(nrSession)])
    calculateCorMats(Params,nrBand,nrSession);
  end
end

% calculate extra statistical maps (t-map, median map, percentile maps):
disp(' ')
disp('Calculating extra statistical maps:')
for nrBand = 0:Pub.nrFreqBands
  for nrSession = 1:Priv.nrSessions
    disp(['Band: ' num2str(nrBand) ', Session: ' num2str(nrSession)])
    calculateStatsMaps(Params,nrBand,nrSession);
  end
end

disp(' ')                                                                                                           
disp(['Calculating inter-subject phase synchronization maps:'])                                                           
for nrBand = 0:Pub.nrFreqBands                                                                                      
  for nrSession = 1:Priv.nrSessions                                                                                 
    disp(['Band: ' num2str(nrBand) ', Session: ' num2str(nrSession)])                                               
    calculatePhaseSynch(Params,nrBand,nrSession);                                                               
  end                                                                                                               
end

% BATCH 4:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% calculate map-wise thresholds using the FDR:
disp(' ')
disp('Calculating critical thresholds:')
for nrSession = 1:Priv.nrSessions
  calculateThresholds(Params,nrSession,0); % across session
  calculateThresholds(Params,nrSession,1); % time-windows
end
% calculate thresholds according to maximal statistic:
calculateThresholdsPF(Params);

% BATCH 5:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Calculate intersubject synchronization curves from time-windowed data:
disp(' ')
disp('Calculating inter-subject synchronization curves:')
for nrBand = 0:Pub.nrFreqBands
  for nrSession = 1:Priv.nrSessions
    disp(['Band: ' num2str(nrBand) ', Session: ' num2str(nrSession)])
    calculateSynchCurves(Params,nrBand,nrSession);
  end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

disp(' ')
disp('Finished!!')
